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Abstract. Meshes composed of well-centered simplices have nice orthogonal dual meshes (the 
dual Voronoi diagram). This is useful for certain numerical algorithms that prefer such primal-dual 
mesh pairs. We prove that well-centered meshes also have optimality properties and relationships 
to Delaunay and minmax angle triangulations. We present an iterative algorithm that seeks to 
transform a given triangulation in two or three dimensions into a well-centered one by minimizing 
a cost function and moving the interior vertices while keeping the mesh connectivity and boundary 
vertices fixed. The cost function is a direct result of a new characterization of well-centeredness 
in arbitrary dimensions that wc present. Ours is the first optimization-based heuristic for well- 
centeredness, and the first one that applies in both two and three dimensions. We show the results of 
applying our algorithm to small and large two-dimensional meshes, some with a complex boundary, 
and obtain a well-centered tetraliedralization of the cube. We also show numerical evidence that 
our algorithm preserves gradation and that it improves the maximum and minimum angles of acute 
triangulations created by the best known previous method. 
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1. Introduction. A completely well- centered mesh is a simplicial mesh in which 
each simplex contains its circumcenter in its interior. A 3-dimensional example is a 
tetrahedral mesh in which the circumcenter of each tetrahedron lies inside it and the 
circumcenter of each triangle face lies inside it. Weaker notions of wcU-ccntcrcdncss 
require that simplices of specific dimensions contain their circumcenters. In two di- 
mensions, a completely well-centered triangulation is the same thing as an acute 
triangulation. 

Typical meshing algorithms do not guarantee well-centeredness. For example, 
a Delaunay triangulation is not necessarily well-centered. In this paper we discuss 
well-centered triangulations, with particular application to triangle and tetrahedral 
meshes. We present an iterative energy minimization approach in which a given 
mesh, after possible preprocessing, may be made well-centered by moving the internal 
vertices while keeping the boundary vertices and connectivity fixed. 

A well-centered (primal) mesh has a corresponding dual mesh assembled from a 
circumcentric subdivision [23] . For an n-dimensional primal mesh, a fc-simplex in the 
primal corresponds to an (n — fc)-cell in the dual. For example, in a well-centered 
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planar triangle mesh, the dual of a primal interior vertex is a convex polygon with 
boundary edges that are orthogonal and dual to primal edges. This orthogonality 
makes it possible to discretize the Hodge star operator of exterior calculus [1] as a 
diagonal matrix, simplifying certain computational methods for solving partial dif- 
ferential equations and for topological calculations. Some numerical methods that 
mention well-centered meshes in this context are the covolume method [29] and Dis- 
crete Exterior Calculus [10, 23]. 

Well-centered meshes are not strictly required for these or other related meth- 
ods; however, some computations would be easier if such meshes were available. For 
example, a stable mixed method for Darcy flow has recently been derived using Dis- 
crete Exterior Calculus [24] and applied to well-centered meshes generated by our 
code and to Delaunay meshes. That numerical method passes patch tests in 2 and 
3 dimensions for both homogeneous and heterogeneous problems. Figure 1.1 (repro- 
duced from [24] by permission of the authors) shows the velocities from a solution to 
the Darcy flow problem in a layered medium. The solution was computed with that 
numerical method and a well-centered mesh. 




Fig. 1.1. Darcy flow in a medium with 5 layers, computed on a well-centered mesh. The odd 
layers have a permeability of 5 and even layers have permeability of 10. The velocities in the odd 
and even layers should be different and should have no vertical component, as shown. The mesh 
was created using our code. Figure taken from [24], used by permission from authors. 



In the case of covolume methods applied to Maxwell's equations, a justification 
for well-centered triangulation is given in [32, 33, 34, 35]. 

Another example from scientific computing is space-time meshing. When tent- 
pitching methods for space-time meshing were first introduced, the initial spatial mesh 
was required to be acute, which for two-dimensional meshes is the same thing as being 
well-centered [38]. More recently this requirement has been avoided, although at the 
expense of some optimality in the construction [20] . 

In two dimensions, well-centered meshes achieve optimality in two objectives that 
are important in some applications. If a planar point set has a well-centered trian- 
gulation, that triangulation both minimizes the maximum angle and maximizes the 
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minimum angle. We don't know any generalizations of this double optimality to 
higher dimensions, but it is known that in any dimension if the convex hull of a point 
set has a well-centered triangulation, then that triangulation is unique and it is the 
Delaunay triangulation [30]. 

2. Our Results. We characterize well-centered triangulations in arbitrary di- 
mensions, prove optimality results for two-dimensional well-centered triangulations, 
and give many experimental results. 

The new characterization of well-centeredness that we give here is a useful the- 
oretical tool that allows us to relate well-centeredness and Delaunay triangulation 
in arbitrary dimensions. In addition, it is also a practical tool since it presents, for 
the first time, a path to the creation of higher-dimensional well-centered triangu- 
lations. Even the formulation of an optimization approach for higher-dimensional 
well-centeredness would be difficult without such a characterization. Indeed, ours 
is the first algorithm to even consider using an optimization approach to seek well- 
centeredness. This approach allows us both to improve existing triangulations in 
and to create well-centered triangulations in MP. We also prove optimality results 
about our cost function and optimality results that relate well-centeredness to well- 
known triangulation schemes. The specific results are enumerated below. 

(a) We introduce a new characterization of well-centeredness in arbitrary dimen- 
sions (Thm. 4.1). (b) As a simple corollary (Cor. 4.2) we show that for any dimension 
n, an n-well-centered triangulation of a convex subset of M" is Delaunay, which is a 
new proof of a result in [30]. (c) Using the characterization of Thm. 4.1 we define a 
family of cost functions Ep (equation 5.2) suitable for creating well-centered triangu- 
lations in arbitrary dimensions, (d) With these we design an algorithm that optimizes 
meshes with the goal of producing well-centered meshes. The algorithm generalizes 
our previous angle-based optimization in two dimensions, described in [40]. Ours is 
the first known strategy for well-centeredness that generalizes to higher dimensions. 

(e) Using the algorithm we produce a well-centered triangulation of a cube (Fig. 7.12). 

(f) We show several two dimensional examples, including one with more than 60000 
triangles (middle of Fig. 7.11). (g) In two dimensions, every algorithm proven to 
generate acute triangulations may produce angles arbitrarily close to 7r/2. More- 
over, in all cases we have tried, our optimization algorithm can improve the quality 
of planar acute-angled triangulations produced by other heuristics for creating acute 
triangulations . A challenging example is shown in Fig. 7.9. (h) We also demonstrate 
numerically that graded triangulations maintain their gradation while being processed 
by our algorithm (Fig. 7.3, 7.8, 7.9). This is useful since producing provably acute 
graded triangulations is an open problem, (i) For planar triangulations, we show 
that the minmax triangulation [17j is the optimal triangulation with respect to our 
energy (Cor. 6.3). (j) We give a different proof for the acute angle case of a result 
from [4]; we show that if a planar point set admits a 2- well-centered triangulation, 
then that triangulation is the unique Delaunay triangulation and the unique minmax 
triangulation of the point set (Thm. 6.4). 

Our experimental results in three dimensions are rudimentary, although even 
these were not available before our work. The difficulty in three dimensions lies fur- 
ther upstream, in a step that precedes the application of our optimization algorithm. 
In the planar case, an interior vertex with four neighbors must be incident to an ob- 
tuse triangle, but some simple connectivity preprocessing can fix this problem [40]. 
Similarly, a tetrahedral mesh may have topological obstructions to well-centeredness. 
The topological obstructions in this case, however, are not yet fully understood. Some 
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progress has been made in our other work [41] by studying the Unk of (topological 
sphere around) a vertex, but much remains to be done. The techniques used to study 
such topological obstructions are interesting, but they are transversal to this paper. 

3. Previous Results. We are concerned with triangulations for which the do- 
main is specified by a polygonal or polyhedral boundary. Our main objective is ob- 
taining well-centered triangle and tetrahedral meshes. Relevant work can be divided 
into constructive and iterative approaches. 

Constructive approaches start with specified input constraints and generate addi- 
tional points, called Steiner points, and/or a corresponding triangulation. Normally 
a point is committed to a position and never moved afterwards. An algorithm for 
nonobtuse planar triangulations based on circle packings is described in [3]. More 
recent works describe improved constructions for nonobtuse triangulations while also 
describing how to derive an acute triangulation from a nonobtuse one [26, 44]. There 
are two major difiiculties with such methods. The first is that these algorithms aim 
to achieve a triangulation of size linear in the input size. As a result, the largest 
and smallest angles can be arbitrarily close to n/2 and respectively. The second 
major diSiculty with these algorithms is that they do not offer a clear path towards 
a higher- dirriensional generalization. Moreover, we are not aware of any existing im- 
plementations of these algorithms, which seem to be primarily of theoretical interest. 
As recently as 2007, Erten and Ungor [21] proposed a variant of the Delaunay refine- 
ment algorithm for generating acute triangulations of planar domains. This heuristic, 
which relocates Steiner points after they are added, has been implemented and ap- 
pears to work quite well. Experiments suggest, however, that the maximum angle in 
the output is often near 7r/2, and our method is able to improve their meshes. See, 
for example, the mesh of Lake Superior in Section 7. 

There is also a constructive algorithm that achieves a well-centered quahty trian- 
gulation of a point set [5] (with no polygonal boundary specified), and an algorithm 
for constructing nonobtuse quality triangulations [27]. Also relevant is an algorithm 
that, given a constraint set of both points and segments in the plane, finds a trian- 
gulation that minimizes the maximum angle [17], without adding points. If an acute 
triangulation exists for the input constraints, the algorithm will find one. The most 
promising of the constructive algorithms is probably [21] mentioned above. But for 
this algorithm, as well as for the others mentioned in this paragraph, we are not aware 
of higher-dimensional generalizations. 

Yet another approach is the mesh stitching approach in [32, 34, 35]. In this 
scheme, the region near the boundary and the interior far from boundary are meshed 
seperately and these two regions are stitched with a special technique. However, in 
three dimensions, the method is unable to generate a well-centered triangulation in 
their examples [32]. 

On the other hand, there are iterative or optimization approaches which allow an 
initial triangulation (possibly the canonical Delaunay) and then move the points while 
possibly changing the connectivity. These algorithms often apply in three dimensions 
as well as two. Moreover, there are many well-known existing meshing algorithms, 
some of which generate quality triangulations [15, 31] and have reliable implementa- 
tions. An iterative approach can start from an existing high-quality mesh and seek 
to make it well-centered whik^ retaining its high quality. 

In the class of iterative approaches there are optimization methods like centroidal 
Voronoi tessellations [12, 13, 14], variational tetrahedral meshing [2]. Each of these 
methods has a global cost function that it attempts to minimize through an iterative 
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procedure that alternates between updating the location of the mesh vertices and the 
triangulation of those vertices. Our algorithm has sonic similarities to these methods, 
but uses a cost function explicitly designed to seek well-centered simplices, in contrast 
to the cost functions optimized in [12] and [2]. 

There are also many iterative optimization methods that, like our mctliocl, relo- 
cate vertices without changing the initial mesh connectivity. Traditional Laplacian 
smoothing [43] is one such method. Such methods improve meshes according to some 
criteria, but do not typically produce well-centered meshes. (See, for example, our 
comparisons with Laplacian smoothing in Sections 7.5 and 7.7.) 

In addition to optimization approaches that work directly with a mesh, there axe 
several algorithms that generate c;irclc pac;kings or circle patterns by optimizing the 
radii of the circles. In particular, the algorithms for creating circle patterns that were 
proposed in [9] and [6] can be adapted to create triangulations. These algorithms 
produce circle patterns that have specified combinatorics, but they do not permit a 
complete specification of the domain boundary, so they are not appropriate to our 
purpose. 

The problem of generating a well-centered tetrahedralization in M.^ is considerably 
harder than the two-dimensional analogue. A complete characterization of the topo- 
logical obstructions to well-centeredness in three dimensions is still an open problem, 
although a start has been made in our work elsewhere [41]. Similarly, the problem 
of generating a three-dimensional acute triangulation — a tetrahedralization in which 
all the dihedral angles are acute — is more difficult than generating a two-dimensional 
acute triangulation. For tetrahedra, it is no longer true that well-centeredness and 
acuteness are equivalent [39, Section 2]. In addition, acute tetrahedralizations are 
known for only restricted domains. For example, until recently it was not known 
whether the cube has an acute triangulation. The construction that showed the cube 
does have an acute triangulation made use of the well-centered optimization discussed 
in this paper [42]. 

4. Characterization of Well-Centeredness. We begin with a new charac- 
terization of well-centeredness in arbitrary dimension. This characterization allows 
us to create an algorithm, described in Section 5, that uses optimization to seek 
well-centeredness. It also serves, later in the current section, as a theoretical tool in 
relating arbitrary-dimensional well-centeredness to Delaunay triangulations. 

Consider an n-dimensional simplex a" embedded in Euclidean space R™, m> n. 
The afhne hull of cr", afF((7"), is the smallest afiine subspace of M"* that contains cr". 
In this case, aff (cr") is a copy of M" embedded in M™. The circumcenter of cr", which 
we denote c(cr"), is the unique point in aff(cr") that is equidistant from every vertex 
of cr". 

For an n-simplex cr" with n > 3, it is possible for cr" to contain its circumcenter 
c(cr") while some proper face -< cr" does not contain its circumcenter c(cr''). It 
is also possible that for all 1 < p < n and all (jP -< cr", c(ct^) lies in the interior of 
(jP, but cr" does not contain its circumcenter. (See [39] for examples with n = 3.) 
Thus we say that an n-simplex cr" is a {pi, ... ,pk) -well-centered simplex if for pi, 
i = 1, . . . , fc, all faces of cr" of dimension Pi < n properly contain their circumcenters. 
The parentheses are suppressed when referring to only one dimension. A simplex cr" 
is completely well- centered if it is (1.2,.... n — 1, n)-well-centered. 

In this section we give an alternate characterization for an n-simplex a" that is n- 
well-centered. The characterization, which shows how the n-well-centered n-simplex 
generalizes the acute triangle to higher dimensions, uses the concept of an equatorial 
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Fig. 4.1. An illustration of the proof of Theorem 4-1 i'>- t^o dimensions. In an n-uiell- centered 
simplex er", vertex Vi and circumcenter c((t") lie in the same open half-space if", the region where 
circumsphere 5"~^ lies outside equatorial ball B" . 

ball, which we now define. 

Let ct" bo a simplex embedded in a hyperplane P"* with m > n. The equatorial 
ball of cr" in P" is the closed ball {x € P" : |x - c(cr")| < i?(cr")}, where c((t") 
is the circumcenter of cr", P((7") its circumradius, and |-| the standard Euclidean 
norm. In this paper we use the notation B((t") for the equatorial ball of cr". The 
notation is used in the context of cr" -< cr"+^, and the hyperplane P™ is understood 
to be afF(cr"+^). The equatorial ball is an extension of the circumball into higher 
dimensions; it is assumed throughout this paper that the circumball and circumsphere 
of a simplex cr" are embedded in aS^(cr"). Note that here and throughout the paper 
we have impHcitly assumed that an n-simplex is fully n-dimensional, though when 
a simplicial mesh is represented on a computer it may be the case that some of the 
simplices are degenerate. 

Theorem 4.1. The n-simplex cr" = voV\ . . . ?;„ is n-well- centered if and only if for 
each « = 0, 1, . . . , n, vertex Vi lies strictly outside P" :~ B{vqVi . . . Vi^iVi^i . . . Vn) ■ 

Proof. Figure 4.1 illustrates this proof in dimension n = 2. It may help the reader 
understand the notation used in the proof and give some intuition for what the proof 
looks like in higher dimensions. 

First we suppose that cr" is n- well-centered. Let 5*"^^ = ^""^(cr") be the cir- 
cumsphere of cr". Now aff(cr") is a copy of M", and within that copy of K", cr" is 
an intersection of half-spaces. Considering some particular vertex of cr", we know 
that one of the bounding hyperplanes of ct" is the hyperplane P"^^ that contains the 
simplex cr""^ = vqVi . . . Vi^iVi+i ...«„. 

Hyperplane P"^^ partitions our copy of M" into two half-spaces — an open half- 
space iJ" that contains the interior of cr" and vertex Vi, and a closed half-space that 
contains cr"~^ (on its boundary). 

Because cr" is well-centered, c(cr") lies in its interior. Thus c(cr") lies in iJ", the 
open half-space that contains u^. Consider, then, the line through c(cr") and c(cr"^^). 
Within if", this line intersects 5"^^ at a point Xi with \xi~c{cr'"') \ — R{a"). Moreover, 
\xi - c(cr"~^)| > P(cr") > i?(cr""^). We see that Xi lies outside Pf and conclude that 
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Fig. 4.2. One characterization of n-well-centeredness of an n-simplex cr^ is that for each 
vertex Vi of a" , Vi lies outside of the equatorial ball B" of the facet it" opposite Vi . 

S""^ n if" lies outside B". In particular, since Vi £ S*""^ fl fff, we know that Vi lies 
outside B". Since Vi was chosen arbitrarily, we conclude that Vi lies outside B" for 
each i = 0, 1, . . . , n, and necessity is proved. 

For sufficiency we consider an n-simplex a" such that Vi lies outside Bf for each 
i = 0, 1, . . . , n. We will show that the circumcenter c((t") lies in the interior of tr" 
by demonstrating that for each vertex Vi, c{a") lies in fff. Wc know that P"~^ cuts 
gn-i jj^^Q pg^j.^ inside B" and a part outside B", and we have just established that 
whichever of the (open) half-spaces contains c(f7") is the half-space where lies 
outside _B". Since wc arc given that Vi G S"~^ lies outside _B", wc know that Vi and 
c(cr") must lie in the same open half-space iJ". This holds for every Vi, so c(c7") is in 
the interior of u", and u" is, by definition, n- well-centered. □ 

Figure 4.2 shows how Thm. 4.1 can be applied to a tetrahedron. In Fig. 4.2 we 
see that for each vertex Vi of the tetrahedron, Vi lies outside of equatorial ball Bf. 
By Thm. 4.1 wc can conclude that the tetrahedron is 3- well-centered, even though we 
have not precisely located its circumcenter. This cleaxly generalizes the acute triangle; 
the angle at vertex Vi of a triangle is acute if and only if Vi lies outside , and a 
triangle is 2-well-ccntcrcd if and only if each of its angles is acute. 

When we say that a mesh is a {pi, . . . ,pk) -well-centered mesh, we mean that every 
element of the mesh is a (pi, . . . ,Pfc)-well-centered simplex. In the proof of Thm. 4.1 
we showed that for each face <j"~^ of an n-wcU-ccntcrcd n- simplex cr". the hypcrplane 
aff(c7"~^) cuts the circumball of cr" into two pieces, one piece contained in B" and 
the other piece lying on the same side of aff(cr"~^) as the interior of ct". It follows 
that the circumball of cr" is contained in (Ui -^?) U cr". (It can be shown, in fact, 
that cr" C Uj Bf , but we do not need that result here.) Moreover, if we consider 
some other n- well-centered n-simplex r" such that cr"~^ = r" fl ct", and if vertex u 
is the vertex of r" opposite (t^~^, then Thm. 4.1 implies that u is outside B". Thus 
u also lies outside the circumball of ct". If the underlying space of the mesh is a 
convex subset of M", we can conclude that the mesh is locally Delaunay. Since in any 
dimension a locally Delaunay mesh is globally Delaunay [16], we obtain a new proof 
of the following result, which was originally proved by Rajan [30]. 

Corollary 4.2. // a simpicial rncsh of a convex subset o/K" is n-well-centered, 
then the mesh is a Delaunay triangulation of its vertices. 

The converse, of course, is not true. Section 6 gives more details for the planar 
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case. 

5. Iterative Energy Minimization. Given a simplicial mesh, we seek to make 
the mesh wcU-cciitcrcd by minimizing a cost function dcfincci over the mesh. We'll 
refer to the cost function as energy. Our method is somewhat similar to the methods of 
[2] and [12] in that it uses an iterative procedure to minimize an energy defined on the 
mesh, but for reasons discussed in Section 6, it differs in that the mesh connectivity 
and boundary vertices remain fixed as the energy is minimized. Also, in contrast to 
the methods of [2] and [12], the cost function we minimize is exphcitly designed to 
achieve the aim of well-centeredness. This section describes the energy we minimize, 
which is the main component of our method. 

Before describing the energy we note that at times the mesh connectivity or 
boundary vertices of an initial mesh are defined in such a way that no well-centered 
mesh exists. For such cases one can apply a preprocessing algorithm to update the 
mesh connectivity. Section 6 discusses this problem in more detail. 

In the proof of Thm. 4. 1 we see that in order for a simplex cr" to be n- well-centered, 
the circumcenter c((7") must lie on the same side of facet cr"~^ as vertex Vi. To convert 
this discrete variable into something quantitative we introduce the function h{vi, cr"), 
the signed distance from c((t") to aff(cr"^^) with the convention that h{vi,a'^) > 
when c((t"') and are on the same side of a,S{a"~^). The magnitude of h{vi,a^) 
can be computed as the distance between c(fT") and c(ct"^^), and its sign can be 
computed by testing whether c(o-") and Vi have the same orientation with respect to 
aff((T"~^). A mesh is n- well-centered if and only if h{vi,a^) > for every vertex Vi 
of every n-simplex cr" of the mesh. 

We divide the quantity h{vi, ct") by the circumradius R{a'^) to get a quantity that 
does not depend on the size of the simplex a". We expect a cost function based on 
h{vi,a"-)/R{cr^) to do a better job than the basic h{vi,a^) at preserving properties 
of the initial mesh. In particular, the grading (relative sizes of the elements) of the 
initial mesh should be preserved better with /i/ J? than with h. Sazonov et al. have also 
noticed that cost functions based on the quantity h/R may be helpful in quantifying 
well-centeredness [32] . 

Note that -1 < /i(i;i,cr")/i?(cr") < 1 for finite ct", because i?(cr")2 = h{v,,a'^f + 
i?(CT"^^)^. Instead of using the quantity h/R directly, we consider the function 



= max 

vertices vGcr" 



Ria^) 



where < fc„ < 1 is a constant that may depend on the dimension n of the simplex. 
The advantage of minimizing /„ as opposed to maximizing h/R is that if fc„ is chosen 

properly, the measure penalizes simplex vertices where h/R approaches 1 (e.g., small 
angles of triangles and sharp points of needle tetrahedra) as well as vertices where 
h/R < 0. 

We want to choose fc„ so that /n(c") is minimized when cr" is the regular n- 
simplex. Taking A:„ = 1/n may seem like a good choice because it is clear that the 
regular simplex minimizes /„. (When fc„ = 1/n, /ra(cr") = for the regular n-simplex 
cr"). We show in Lemma 5.1, however, that the regular simplex minimizes /„ for any 
l>kn> 1/n. 

Lemma 5.1. For kn > 1/n, the measure fnicr"") is minimized when a" is a regular 
simplex. 
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Proof. Suppose that kn > For the regular simplex, then, /ri(cr") = kn — 1/n. 
Thus it suffices to show that for any simplex tr" there exists a vertex v such that 
h{v,a"') < R{a")/n; at such a vertex we have 



h{v,a") 



i?((T") 



Wc have seen that for a simplex that is not n-well-centered, there exists a vertex 
V with h{v, tr") < 0, so it remains to prove this for simplices that are n-wcU-centered. 

Suppose a" is n-well-centered. Let h := miuj h{vi, ct"). Consider a sphere S^~^ C 
aff(cr") with center c((t") and radius h. Wc claim that ct" contains the sphere S"""^. 
Indeed, for each facet <j"^^ of cr", since the radius of S*""^ is ft- < h{vi,a") we have 
that the sphere S'^~^ is contained in the same half space as c((t") and Vi. Thus the 
sph(^r(^ is contained in the intersection of half spaces that defines the simplex, i.e., is 
contained in the simplex. 

It follows, then, that h < r{a") where r{a") is the inradius of ct". We know that 
h/R < r/R < 1/n and that equality is achieved for only the regular simplex. (The 
inequality r/R < 1/n is proved in [25], among others.) □ 



In light of Lemma 5.1, taking fc„ = 1/2, independent of n, is a good strategy, 
because for kn = 1/2 the cost function will prefer any n-well-centered simplex to 
any simplex that is not n-well-centered, and among all n-well-centered simplices, /„ 
will prefer the regular simplex over all others. We use fc„ = 1/2 for all of the results 
discussed in Section 7. 

For kn > the objective of n-well-centeredness is achieved when \h/R — kn\ < kn 
at every vertex of every simplex cr". (Note that this is not a necessary condition 
if kn < 1/2.) Our goal, then, is to minimize \h/R—kn\ over all vertices and all 
simplices, driving it below fc„ at every vertex of every simplex. It could be effective 
to work directly with 



(V,T)= max 

simplices ct^GT 
vertices Vi^a'^DV 



h{vua^) _ 1 

i?(CT") 2 



(5.1) 



but we choose instead to minimize an approximation to 2£oo given by 



E^{M) = E^{V,T)= J2 



2h{vi,a") 



Ria"") 



- 1 



(5.2) 



where p is a parameter. A4 here stands for a mesh consisting of vertices V with 
particular coordinates and a connectivity table T that describes which groups of 
vertices form simplices. Note that linip^oc {Ep (M))^^^ = 2Eoo{M), so Ep{M.) is 
indeed an approximation to 2Eo^{M.). The factor of 2 is included for numerical 
robustness. The parameter p influences the relative importance of the worst vertex- 
simplex pair compared to the other vertex-simplex pairs in computing the quality of 
the mesh as a whole. It is convenient to choose p as a positive even integer, since the 
absolute value need not be taken explicitly in those cases. 

As stated, the measure Ep{M.) leaves some ambiguity in the case of a degenerate 
simplex, which may occur in a computational setting. For several reasons, including 
a desire to maintain upper semicontinuity of the cost function, we use the convention 
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Fig. 5.1. For a triangle, h/R = cos{9). 



that any degenerate simplex, even one with coincident vertices, has its circumcenter 
at infinity and h/R = —1. 

Figure 5.1 shows the quantities h and ii in a sample triangle. We see in the figure 
that cos(^) = h/R. Thus (5.2) is a generalization of the energy 



which is a constant multiple of the energy the authors proposed earlier for achieving 
well-centeredness of planar triangle meshes [40]. In three dimensions the quantity 
h/R is related to the cosine of the tetrahedron vertex angle, as discussed in [32]. 

The cost functions Ep and E.^ are not convex. When designing a cost function 
for mesh optimization, one might hope to develop a function that is convex, or, if 
not convex, at least one that has a unique minimum. It is, however, not possible to 
define an energy that accurately reflects the goals of wcU-ccntcred meshing and also 
has a unique minimum. Consider the mesh shown on the left in Fig. 5.2. We suppose 
that the boundary vertices are fixed, but the interior vertex is free to move. We want 
to decide where to move the interior vertex in order to obtain a well-centered mesh. 
The right side of Fig. 5.2 shows where the free vertex can be placed to produce a 
well-centered mesh. The light gray regions are not allowed because placing the free 
vertex in those regions would make some boundary angle nonacute. (The dotted lines 
indicate how the four most important boundary angles infiuence the definition of this 
region.) The darker gray regions, shown overlaying the light gray region, are not 
permitted because placing the interior vertex in those regions would make some angle 
at the interior vertex nonacute. 

If the interior vertex is placed in either of the two small white regions that remain, 
the mesh will be well-centered. We see that the points permitted for well-centeredness 
form a disconnected set in M^. Moreover, the mesh is radially symmetric, so there 
is no way to create an energy that prefers one white region over the other unless we 
violate the desired property that the energy be insensitive to a rotation of the entire 
mesh. Any symmetric energy that has minima in only the white regions must have 
at least two distinct global minima. 

In most planar triangle meshes there is an interior vertex v that has exactly six 
neighbors, all of which are interior vertices. If all interior vertices are free to move, 
as we assume in the method we propose, then the six neighbors could be moved 
into the relative positions that the boundary vertices have in the mesh in Fig. 5.2. 




(5.3) 



eeM 
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Fic;. 5.2. A cost function that accurately reflects the goal of well-centeredness cannot have a 
unique minimum, because the set of points that make the mesh well-centered may be a symmetric 
disconnected set. 



Moving V around when its neighbors have such positions should exhibit nonconvexity 
in whatever cost function we might define. 

6. The Optimal Planar Triangulation. A variety of our experimental results 
appears in Section 7 below. The results support the claim that Ep is an appropriate 
cost function for quantifying the 2-well-centeredness of a planar mesh. In some cases, 

though, the mesh connectivity, the fixed boundary vertices, or a combination of the 
two are specified in such a way that no well-centered mesh exists with the given mesh 
connectivity and boundary vertices. The simplest example of this is a planar mesh 
with an interior vertex v that has fewer than five neighbors. Since the angles around 
V sum to 27r, v has some adjacent angle of at least 7r/2. The triangle containing 
that angle is not 2- well-centered. Similarly, a boundary vertex with a boundary angle 
measuring at least tt/2 must have enough interior neighbors to divide the boundary 
angle into pieces strictly smaller than w/2. We will refer to a vertex that does not 
have enough neighbors as a lonely vertex. (In three dimensions, a vertex must have 
at least 7 incident edges to permit a 3-well-centered mesh, though having 7 neighbors 
is not sufficient to guarantee that a 3-well-centered neighborhood exists.) 

One way to approach problems with mesh connectivity, such as the problem of 
lonely vertices, is a global mesh connectivity update, i.e., to change the mesh con- 
nectivity over the entire mesh. The methods that use Voronoi diagrams [12] and 
variational triangulations [2] both employ this approach, updating to a Delaunay 
mesh each time the vertices are relocated. In this section we show that the optimal 
triangulation of a planar point set with respect to the energy E^c is a minmax trian- 
gulation, i.e. a triangulation that minimizes the maximum angle. Note that in general 
a minmax triangulation is not a DelaTinay triangulation. (A Delaunay triangulation 
is, rather, a maxmin triangulation of a planar point set [37]). 

There is an 0{n^ logn) time algorithm for computing the minmax angle triangu- 
lation of a fixed set of points in the plane [17], so in the plane it might be feasible to 
recompute the optimal triangulation at every step of our iterative algorithm. It is not 
clear, however, whether the algorithm of [17] can be generalized into higher dimen- 
sions. At the end of this section we discuss some other reasons to avoid recomputing 
the optimal triangulation after each step of energy minimization. 

In the rest of this section we restrict our attention to a given set of vertices V 
in M^, fixed at their initial locations. Given V we seek the mesh connectivity T that 
minimizes Eao(y,T). Throughout this section, where we refer to mesh connectivity 
or triangulation it is assumed (often implicitly) that we mean an admissible triangu- 
lation, i.e., a triangulation of V that covers the convex hull of V, conv(V), and has no 
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inverted or overlapping triangles. Many of the results would apply when considering 
a different set of admissible triangulations, but some might need small modifications, 
depending on the particular set of triangulations admitted. 

Since we are working in the plane, the discussion is based on planar angles 6 and 
the cost function defined in (5.3) in terms of cos(^). In particular we consider the cost 
functions 



where in the latter two cases we require 9 G [0,n]. 

We start by showing that when all triangulations of a planar point set have 
a maximum angle that is at least tt/2, a triangulation minimizing Emax is also a 
triangulation that minimizes -Ecos- This claim is readily proved as a corollary of the 
following proposition. 

Proposition 6.1. Let f be a strictly increasing function of 6 and g a nondecreas- 
ing function of 6 for G [0,7r]. If Ef{T) = max{/(6li)} and Eg{T) = max{5(6»i)}, 
then arg min Ef C arg min Eg . 

Proof. For each triangulation T, there exists some angle Or such that Ef{T) = 
max{/((?i)} = f{Or)- Thus for all other angles 6 appearing in triangulation T, we 
have that /(6»r) > f{0)- 

Consider a specific triangulation % e arg mini?/. We have Ef{To) < Ef{T) for 
all triangulations T. Thus /(^Tq) < /(^r) Moreover, since / is a strictly increasing 
function of 6, we can conclude that Otq < Or Then since g is nondecreasing, we have 
gi^To) < 9{9t) for all triangulations T. 

Now we claim that for arbitrary triangulation T we have g{Or) > g{6) for all 
angles 9 appearing in triangulation T. If this were not the case, then there would exist 
some angle ^ in T with g{9) > g{9r)- Since g is nondecreasing, it would follow that 
6 > 9t, and since / is strictly increasing, we would have f{9) > f{9r)- This, however, 
contradicts our definition of 9r, which states that /(^t) = iiiax{/(^i)} > f(9). We 
conclude that the claim is correct. 

It follows, then, that g{9r) = max{5f(^j)} = Eg(T) for each triangulation T. 
In particular, the inequality g{9ro) < 9{9t) implies that Eg{To) < Eg{T) for all 
triangulations T. By definition, is a. member of the set arg min E^^. □ 

Corollary 6.2. If f is a strictly increasing function of 9 for 9 e [0,7r], then 
arg min Ej = arg min Emax ■ 

Proof. The function E^ax is of the form Eg where g is the identity function 
on [0,7r]. Since g is a strictly increasing function, we may apply Proposition 6.1 in 
both directions to show that axginin Ef C aigmin E^ax and that argmin E^ax Q 
axgrnin Ef. We conclude that a,j:gmm Emax = argminiJ/. □ 

Corollary 6.3. // all triangulations of a set of vertices V that cover conv(V) 
have maximum angle at least 7r/2, then a triangulation minimizing E^ax also mini- 
mizes Ecos o-nd vice versa. 




Emin (V, T) 




Emax (V, T) 
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Proof. We can restate the corollary as follows. If Emax > 7r/2 for all triangulations 
T, then argminii'cos = arg min Sm^K • This follows because Ecos is of the form Ef 
where / = |2cos(0) — 1| is a strictly increasing function on the interval [7r/2,7r], and 
f{0) < /(7r/2) for < ^ < 7r/2. For all practical purposes, we could redefine / on 
[0, 7r/2) to make / a strictly increasing function on [0, tt]. The redefinition would have 
no effect because for all T, the maximal f{9i) occurs at some 9i > tt/2. 

Some care should be taken if we allow meshes that have an angle ^ = 0, but we 
know that a triangle with an angle of has some angle measuring at least tt/2, even 
if two of the triangle vertices coincide. Since /(7r/2) = /(O), we may say that on a 
triangle with angle 0, / is maximized at the largest angle 6 > 7r/2. □ 

It should be clear that the proofs of Prop. 6.1 and Cor. 6.3 do not apply when a 
triangulation exists with Emax < 7i'/2. In that case, Ecos may be maximized at some 
angle 9^0 rather than at the largest angle of the mesh. In the next theorem we es- 
tablish that there is an important relationship between axg min E^ax and arg min E^-^g 
even when a well-centered triangulation exists. (This theorem is the acute angle case 
of a result from [4], presented here with a different proof.) 

Theorem 6.4. If a 2-well- centered triangulation of a planar point set exists, 
then that 2-'well- centered triangulation is unique and is both the unique Delaunay 
triangulation of the point set and the unique minmax triangulation of the point set. 

Proof. Recall that if the Delaunay complex of a planar point set has a cell that is 
not triangular, then this cell is a convex polygon with more than three vertices. The 
vertices of the polygon axe all cocircular, and the circumcircle is empty of other points. 
In this case a (nonunique) Delaunay triangulation may be obtained by triangulating 
each such polygon arbitrarily. Any such Delaunay triangulation must contain an angle 
with measure 7r/2 or larger. 

This can be argued from considering the possible triangulations of a Delaunay cell 
that is not triangular. An ear of the triangulation of the Delaunay cell is a triangle 
bounded by one diagonal and two edges of the Delaunay cell. Since the Delaunay cell 
has four or more vertices, at least two triangles will be ears in any triangulation of 
the cell. Moreover, we can divide the circumdisk of the Delaunay cell into a pair of 
closed semidisks in such a way that at least one semidisk completely contains an ear. 
In an ear contained in a semidisk, the angle along the boundary of the Delaunay cell 
is at least tt/2. We conclude that if the Delaunay complex of a planar point set is not 
a triangulation, then no completion of the Delaunay complex to a triangulation (i.e., 
a Delaunay triangulation) yields a 2-well-centered triangulation. 

Suppose, then, that a point set permits a 2-well-centered triangulation Tq. By 
Cor. 4.2, 7q is a Delaunay triangulation. The Delaunay triangulation is unique in 
this case (by the argument of the preceding paragraph). Moreover, any other trian- 
gulation T of the point set has a maximum angle that is at least as large as tt/2. 
(If not, T would be 2-well-centered, and, therefore, a Delaunay triangulation, contra- 
dicting the uniqueness of the Delaunay triangulation.) We conclude that the minmax 
triangulation in this case is Tq and is unique. □ 

Combining Thm. 6.4 with Cor. 6.3 we see that arg min Ef^Qg — arg min E^ax m 
all cases. 

Unfortunately, the minmax triangulation and the Delaunay triangulation both 
have the undesirable property that they may have interior vertices with only four 
neighbors, i.e., lonely vertices. Figure 6.1 shows a small point set for which the 
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Fig. 6.1. The minmax triangulation may produce a triangulation in which interior vertices 
are lonely, even when there are triangulations with no lonely vertices. The sequence of figures 
shows a point set, the minmax triangulation of the point set, an alternate triangulation of the point 
set with no lonely vertices, and a 2-well-centered triangulation that is obtained from the alternate 
triangulation by optimizing E4. 



minmax triangulation contains an interior vertex with only four neighbors. In this 
particular case, the minmax triangulation gives a mesh for which the vertex locations 
optimize both E^. and E/^. Thus optimizing i?oo or E^ will not change the mesh, even 
if we interleave the mesh optimization with recomputing the optimal triangulation. 

As long as we maintain the mesh connectivity given by this minmax triangulation, 
we cannot make the mesh 2-wcll-ccntcrcd, regardless of what function wc optimize. 
To address this problem we choose to use an algorithm that preprocesses the mesh, 
updating the mesh connectivity locally to eliminate lonely vertices. The algorithm we 
use for the two-dimensional case is outlined in [40] . The preprocessing step applied to 
the minmax triangulation produces an alternate triangulation of the initial vertex set. 
(See Fig. 6.1.) For the new triangulation, optimizing E4 quickly finds a 2-well-centered 
mesh. 

A key reason that we choose to preserve the mesh connectivity throughout the 
optimization process is that we want to prevent the appearance of lonely vertices 
during the optimization process. It might be interesting to interleave the energy 
optimization with a retriangulation step that computed a triangulation that minimizes 
the maximum angle among all triangulations with no lonely vertices, but we do not 
know how to compute such a triangulation efficiently. The choice to maintain mesh 
connectivity during optimization also simplifies the handling of meshes of domains 
with holes. 

7. Experimental Results. In this section we give some experimental results 
of applying our energy minimization to a variety of meshes. All of the initial meshes 
shown here permit well-centered triangulations, in many cases because the "initial 
mesh" is the output of some preprocessing algorithm that improves the mesh con- 
nectivity, e.g., the preprocessing algorithm described in [40]. The mesh optimization 
was implemented using the Mesquite library developed at Sandia National Labora- 
tories [8]. We implemented the cost function Ep by writing a new element-based 
QualityMetric with a constructor ac;c;epting the argument p and summing the en- 
ergy values on each element with the standard LPtoPTemplate objective function 
(with power 1). 

We used Mesquite's implementation of the conjugate gradient method to optimize 
Ep on each mesh shown. We did not write code for an analytical gradient, so Mesquite 
numerically estimated the gradients needed for the conjugate gradient optimization. 
The optimization was terminated with a TerminationCriterion based on the number 
of iterations, so where the phrase number of iterations appears in the experimental 
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Fig. 7.1. For two-dimensional meshes, the shade of a triangle indicates the measure of its 
largest angle. 




Fig. 7.2. Prom the initial mesh shown at left, with 3.10% of its triangles nonacute, minimizing 
Ei produces the 2-well-centered mesh shown at right in 30 iterations. Histograms of the angles in 
the mesh are included, with the minimum and maximum angles marked on each histogram. The 
optimization took 1.61 seconds. 

results, it refers to the number of iterations of the conjugate gradient method. For 

the three-dimensional meshes shown here we used the cost function Ep for dimension 
n = 3, which is designed to find 3-well-centered meshes and is not sensitive to whether 
the facets of the tetrahedra are acute triangles. 

All of the experimental results discussed in this section were run on a desktop 
machine with a dual 1.42 GHz PowerPC G4 processor and 2 GB of memory. As is 
often the case with mesh optimization, the algorithm is quite slow. There are cer- 
tainly opportunities for improving the efficiency of the algorithm as well; the authors 
suspect that modifying the algorithm to do optimization only in the regions where it 
is necessary, instead of optimizing over the entire mesh, could improve the efficiency 
significantly. 

Shading scheme: For all the two-dimensional meshes shown in this section, we 
use the scale shown in Fig. 7.1 to determine the shade of each triangle. The shade 
of a triangle is determined by the measure of the largest angle of the triangle. The 
shade gets darker as the largest angle increases, with a noticeable jump at 90° so that 
2-well-centered triangles can be distinguished from nonacute triangles. For example, 
the three meshes in Fig. 6.1 use this shading scheme, and it should be easy to identify 
the triangles that are not 2-well-centered in the first two meshes. 

Along with figures of meshes, we include histograms that show the distribution of 
the angles for two-dimensional meshes. We report near the histogram the percentage 
p and the number n of nonacute triangles in each mesh. The mean of each distribution 
is 60°, and the standard deviation a is written near the distribution. 

7.1. Mesh of a Disk. The mesh of the disk in Fig. 7.2 is small enough that the 
results of an experiment on the mesh can be visually inspected. Many of the triangles 
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Fig. 7.3. Results of an experiment with a mesh of a 2-dimensional slice of the combustion 
chamber inside the Titan IV rocket. The initial rnesh is displayed at the top. Below it is the result 
mesh, which was obtained by 1000 iterations minimizing Eio on the mesh. Histograms show the 
distribution of angles in the initial and final meshes. The zoomed in views of the joint slot ( at the 
top center of the full mesh) show the level of mesh refinement in the regions of higher detail. For 
the histograms and the zoom.ed views, the original m,esh is on the left, and the result mesh is on the 
right. The optimization took 805.35 seconds. 

are already acute in the initial mesh, but some are not. Based on the shading scheme, 
we see visually that the result mesh has no nonacute triangles. The histograms of the 
angles in the mesh confirm this, showing that the maximum angle was reduced from 
121.22° to 82.55°, and the minimum angle has increased from 22.15° to 33.46°. The 
optimization took 1.61 seconds. 

7.2. A Larger Mesh. In Fig. 7.3 wc show results for a larger mesh, a mesh of 
a two-dimensional slice of the combustion chamber inside the Titan IV rocket. This 
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mesh, which is based on a mesh that the third author produced from his work for 
the Center for Simulation of Advanced Rockets, has 8966 triangles. At the top of 
Fig. 7.3 we show an overview of the entire mesh, with the initial mesh at the very top 
and the result (after optimizing Eio for 1000 iterations) just below it. These meshes 
are drawn without showing element edges, because even the thinnest possible edges 
would entirely obscure some parts of the mesh. The background color helps define 
the boundary of the mesh by providing contrast with the light gray elements. 

Below the mesh overview is a zoomed view of the top center portion of the mesh, 
which represents a portion of a joint slot of the titan IV rocket. Figure 7.3 also includes 
histograms of the angle distribution of the full mesh before and after the optimization. 
The angle histogram and zoomed portion for the initial mesh are shown on the left, 
and for the optimized mesh are shown on the right. 

In the initial mesh there are 1188 nonacute triangles (« 13.25% of the triangles), 
with a maximum angle around 155.89°. The result mesh has a maximum angle of 
89.98°, and all but 143 triangles (« 1.59%) have maximum angle below 85°. Of the 
143 triangles that have angles above 85°, 14 have all three vertices on the boundary 
and are thus completely specified by the boundary. One example of this is in the 
upper left corner of the zoomed view, where there is a triangle that looks much like 
an isosceles right triangle. Another 60 triangles are forced to have triangles larger than 
85° because they are part of a pair of triangles along a part of the boundary with small 
but nonzero curvature. There are four such pairs along each curved boundary in the 
zoomed view in Fig. 7.3. In fact, all but 4 of the 143 "worst" triangles have at least 
one boundary vertex, and the remaining 4 triangles each have a vertex that is distance 
one from the boundary. 

7.3. Some More Difficult Tests. The next mesh is a mesh of a circular domain 
with two circular holes. The initial mesh is far from being 2-well-centered, with 
61.04% of its triangles nonacute, and a standard deviation a « 31.238 for the angle 
distribution. An initial attempt to make the mesh well-centered was unsuccessful, but 
two slightly diS^erent strategies, described later, do produce a well-centered mesh. The 
initial mesh and its angle histogram are shown in Fig. 7.4 (left) along with the result 
of minimizing E4 on the mesh for 500 iterations (right). In this case, the optimization 
took 88.70 seconds. Comparing the optimized mesh to the initial mesh we see that the 
quality has improved; the percentage of nonacnite triangles is reduced, the standard 
deviation has improved, and many of the largest angles have been reduced. 

Unfortunately, some of the smallest angles of the initial mesh have also gotten 
smaller in the optimized mesh. In fact, four angles got so small that their triangles 
became inverted in the optimized mesh. The inverted triangles are too thin to actually 
see, but there is one pair near the top right of the mesh and one pair near the bottom 
left. The energy value required to invert a triangle is fairly large, but for large meshes 
or meshes with a high percentage of bad triangles, improvements at other locations 
in the mesh may be significant enough to overcome the cost of triangle inversion for 
a small number of the triangles in the mesh, and using the basic energy Ep can lead 
to inverted triangles. Triangle inversion can be prevented by including an inversion 
barrier in the cost function. 

Energy combined with inversion barrier. Modifying the energy by introduc- 
ing a term that has a barrier against inversion, i.e., a term for which the energy value 
goes to infinity as a triangle moves towards becoming degenerate, is probably the best 
way to handle the problem of triangles that would become inverted with the basic 
Ep. The IdealWeightlnverseMeanRatio QualityMetric provided by Mesquite is a 
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Fig. 7.4. A first attempy at energy minimization applied to the two holes mesh on the left 

does not yield a well-centered mesh. Result after 500 iterations of E4 minimization is shown on the 
right. The optimization took 88.70 seconds. The result mesh has some inverted triangles which are 
too thin to be seen. In subsequent figures we show several strategies for producing a well-centered 
configuration. 




Fig. 7.5. A 2-well-centered mesh of the two holes domain conforming to the mesh connectivity 
and boundary vertices of the original two holes mesh shown in Fig. 7.4- The mesh was obtained using 
slightly modified cost functions Ep that have a barrier against triangle inversion. The optimization 
procedure was 500 iterations of E4, followed by 500 iterations of Eg followed by 500 iterations of Eio ■ 
Total optimization time was 115.37 seconds. 

cost function that has an imphcit barrier against inversion [28]. Let -Eimr represent 
the cost function associated with the IdealWeightlnverseMeanRatio. One can take 
a linear combination of the energy Ep with Eimr to create a new energy that has a 
barrier against inversion and, depending on the coefficients, is still very much like Ep. 
We have found that the energy Ep := lOOi^p + i^imr is often effective in cases where 
the basic Ep leads to inverted triangles. For this problem, for example, using Ep gives 
a well-centered result with no inverted triangles. Starting from the initial mesh and 
applying 500 iterations of E4 followed by 500 iterations of Eq and 500 iterations of £10 
produced the 2-well-centered mesh of the original domain displayed in Fig. 7.5. The 
optimization took 37.37 + 36.79 + 41.21 = 115.37 seconds. 

Improved boundciry vertex locations. Another way to get a well-centered 
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Fig. 7.6. This mesh has the same mesh connectivity as the initial mesh in Fig. 7-4, but the 
vertices along the boundary (and in the interior) have been moved. The 2-well-centered mesh on the 
right was obtained in 18.03 seconds with 200 iterations of Eg minimization. 

mesh from this initial mcsli is to make the optimization problem easier by changing 
the location of the boundary vertices. The mesh on the left in Fig. 7.6 has the same 
mesh connectivity as the initial two holes mesh from Fig. 7.4, but the vertices along 
the boundary have moved. In the initial mesh the vertices along each boundary were 
equally spaced, but in this case, the vertices on the outer boundary are more dense 
at the north and south and less dense at the east and west. The vertices along the 
inner boundary curves have also moved a bit. For this mesh we use the basic energy 
Eq, reaching a well-centered configuration by 200 iterations. The result, obtained in 
18.03 seconds, appears on the right in Fig. 7.6. 

Different mesh of the same domain. The difficulty of finding a 2-well- 
centered mesh is primarily due to the combined constraints of the mesh connectivity of 
the initial mesh and the locations of the boundary vertices. The shape of the domain 
or the fact that the domain is not simply connected axe not inherently difficult for the 
problem of 2-well-centered triangulation. When separated from the mesh connectivity 
of the initial mesh, the location of the boundary vertices are not a problem either. We 
demonstrate this by an experiment on the same domain with a completely different 
mesh that has the same set of boundary vertices and the same boundary vertex loca- 
tions as the meshes of Figs. 7.4 and 7.5. The experiment, shown in Fig. 7.7, produced 
a mesh of the domain with maximum angle around 79.50° by optimizing Eg for 100 
iterations. The optimization took 7.44 seconds. 

7.4. A Graded Mesh. The two holes mesh of Fig. 7.4 and the mesh in Fig. 7.3 
related to the titan rocket are both graded meshes. However, the gradation of those 
meshes was controlled partly by the size of elements on the boundary and by the 
geometry of the mesh. In Fig. 7.8 we show the results of applying energy minimization 
to a mesh of the squaxe with an artificially inducied gradation. The initial mesh and 
angle histogram appear at left in Fig. 7.8. The nearly converged result produced by 
30 iterations minimizing E4 is displayed to its right. 

The initial size of the triangles of a mesh is not always preserved well when 
optimizing the energy. We expect, however, that the energy will generally preserve 
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Fig. 7.7. This is a mesh with the same domain and same boundary vertices as the mesh in 
Fig. 7.4- The 2-well-centered mesh on the right was obtained from the mesh on the left in 7.44 
seconds by minimizing Eg for 100 iterations. The high-quality result shows that the difficulty of 
getting a 2-well-centered mesh in Fig. 7.4 is not due solely to the domain or the boundary vertices. 
The initial mesh for this experiment was generated using the freely available software Triangle [36] 
and heuristics for improving the mesh connectivity. 

the grading of an input mesh if the initial mesh is relatively high quality. This 
hypothesis stems from the observation that the energy is independent of triangle size, 
the idea that the mesh connectivity combined with the property of 2-well-centeredness 
somehow controls the triangle size, and the supporting evidence of this particular 
experiment. 

Thus optimizing graded meshes is a useful application of our algorithm; there are 

no known provably correct algorithms for creation of graded acute-angled triangu- 
lations of planar domains. The recent algorithm of [21] has produced graded acute 
triangulations in a variety of experiments, but in all cases we have tried, we have 
been able to improve the quality of their triangulations (Section 7.5). Moreover, their 
algorithm is not known to generalize to higher dimensions. 

7.5. Mesh of Lake Superior. The Lake Superior domain, with its complicated 
shape, has appeared in many papers about quality meshing. We include an example 
optimizing a mesh of this well-known cioniain. The initial mesh is already 2-well- 
centered in this experiment, but we show that we can improve its quality with our 
optimization algorithm. The results are represented graphically in Fig. 7.9. 

The initial acute-angled mesh is from the work of Erten and Ungor [21] on gen- 
erating acute 2-D triangulations with a variant of Delaunay refinement. The initial 
mesh has a maximum angle of 89.00° with 174 triangles having angles larger than 
88.00°. Directly optimizing £^10 on the initial mesh, Mesquite finds a local minimum 
of Eio after 6.63 seconds (21 iterations). The local minimum has exactly one nona- 
cute triangle (maximum angle 91.03°) and only 40 triangles having angles larger than 
88.00°. The angle histogram for this result is included in Fig. 7.9 at top center. The 
mesh is visually very similar to th(^ initial mesh and does not appear in this paper. 

If we start by optimizing E4 and follow that by optimizing Eiq we obtain a local 
(perhaps also global) minimum of Eiq with with much lower energy than the result 
obtained by directly optimizing Eiq. The result of this optimization process is shown 
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Fig. 7.8. For this graded mesh of the square, minimizing E4 on the initial mesh (left) produces 
a 2-well-centered mesh (right) that has grading similar to the initial mesh. The optimization ran 
for 30 iterations, completing in 2.16 seconds. 

on the right in Fig. 7.9. The optimization took 131.48 seconds total; Mesquite spent 
102.81 seconds (453 iterations) finding a minimum of E4 and 28.67 seconds (125 
iterations) finding a minimum of £^10. 

Laplacian smoothing is a popular mesh optimization technique that was first used 
for structured meshes with quadrilateral elements and later generalized to triangle 
meshes [43]. A brief description of Laplacian smoothing is given in [22]. We compare 
our mesh optimization technique with Laplacian smoothing, using the implementation 
of Laplacian smoothing provided by the Mesquite library. The result of Laplacian 
smoothing on the Lake Superior mesh is shown in Fig. 7.10. The optimization was 
terminated after 100 iterations, which is near convergence. The run time was 1.31 
seconds. The maximum angle in the result is 109.27° and more than 4% of the triangles 
arc nonacutc. 

The result of optimizing the Lake Superior mesh with Laplacian smoothing is 
typical of the results obtained with Laplacian smoothing. We performed experiments 
with Laplacian smoothing on all of the 2-D meshes presented in this paper, and 
no mesh became well-centered except for the mesh of the square in Fig. 7.8, where 
Laplacian smoothing produced a mesh with maximum angle 87.54° compared to the 
max;imum angle of 78.50° obtained by our method. In most cases the percentage 
of nonacute triangles after Laplacian smoothing was between 1% and 5%, but for 
the meshes in Figs. 7.4, 7.5, and 7.6, the percentage of nonacute triangles was much 
higher, getting as high as 48.70% for the mesh in Fig. 7.6. Clearly the traditional 
Laplacian smoothing is not an appropriate tool for finding acute triangulations. 

7.6. Colombia, India, and Thailand. We end om- 2-D experimental results 
with a collection of three large meshes of complicated geographical domains. The 
experiments are summarized in 7.11. For each of these meshes the optimization started 
by minimizing Es for 500 iterations and then proceeded by minimizing Eg, running 
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Fig. 7.9. Result for a mesh of Lake Superior. The initial mesh shown on the left is a 2-weU- 
centered mesh from [21]. The improved mesh shown on the right was obtained by first optimizing 
E4 and then optimizing Eiq. The angle histogram at top center shows the result of optimizing Eio 
directly on the initial mesh. Many of the angles that were near 90 ° have dropped to below 80 °. 




Fig. 7.10. Result of applying Laplacian smoothing to the initial acute mesh of Lake Superior 
(left side of Fig. 7.9). More than 4% of the triangles become nonacute, and the maximum angle 
increases to 109.27°. 

500 iterations at a time until the mesh became well-centered. After the mesh became 
well-centered, we used one more round of 500 iterations minimizing Eg. to get some 
additional improvement in the angle distribution. 

The total number of iterations for the meshes was 2000 iterations for Colombia, 
3500 iterations for India, and 3000 iterations for Thailand, with total optimization 
times of 5284.01 seconds, 16162.20 seconds, and 8263.82 seconds. The meshes are 
quite large, with 38233 triangles, 62370 triangles, and 34562 triangles respectively. In 
each case, more than 19% of the triangles are nonacute in the initial mesh, and the 
maximum angle is larger than 160°, yet the optimization finds a well-centered result. 
It is also clear that the optimization preserves the gradual change in element size from 
the tiny triangles needed to resolve the boundaries to the much larger triangles in the 
interiors of the meshes. 

7.7. 3D Meshes. For tetrahedral meshes, the question of whether the mesh 
connectivity permits a well-centered mesh is more difficult than its two-dimensional 
analogue [41]. In part because we do not yet have an effective preprocessing algorithm 
for tetrahedral meshes, many of our optimization experiments in three dimensions 
have been limited to meshes with carefully designed mesh connectivity. The mesh 
shown in Fig. 7.12 is one of these meshes. The shading of the tetrahedral elements 
in Fig. 7.12 represents the shadows that would result from viewing the faceted object 
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Fig. 7.11. These meshes of complicated geographical boundaries were optimized with an initial 
500 iterations of Eg, followed by successive rounds of 500 iterations minimizing the basic Eg. The 
optimization produces well-centered meshes that preserve the grading of the input meshes. (The dark 
regions near the boundaries of the meshes come from the agglomeration of the edges of triangles 
that are too small to be seen.) The data for the geographical boundaries was produced using the 
CountryDataf] command of Mathematica. Initial meshes were constructed from the input polygons 
using Triangle [36] and heuristics for improving the mesh connectivity. 



under a light source; it has nothing to do with the quahty of the elements of the 
mesh. The full mesh is a mesh of the three-dimensional cube with 430 tetrahedra. 
Figure 7.12 uses a cutaway view to display some of the elements in the interior of the 
mesh. 

Although the initial mesh was carefully designed to have good mesh connectivity 
(e.g., each vertex has at least 10 incident edges) and a high-quality surface mesh, it 
was not 3- well-centered. In fact, 22.33% of the tetrahedra are not 3- well-centered. Op- 
timizing EiQ for 3.92 seconds (20 iterations) produced a 3-well-centered mesh. Even 
though the initial mesh was carefully designed, the optimization result is nontriv- 
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Fig. 7.12. A cutout view showing the interior of a ^-well-centered mesh of the cube. The 
mesh is the result of 3.92 seconds (20 iterations) of optimizing Eiq on an initial mesh for 22.33% 
of the tetrahedra were not Z -well- centered. Recall that a tetrahedron cr^ is 3-well-centered if and 
only if h{v, tJ^)/ R{(^^) > for each vertex v of <t^ . For a regular tetrahedron, h/R= 1/3. The h/R 
distributions for the initial mesh, the result of optimizing Eig , and the result of Laplacian smoothing 
show the superiority of our method for finding 3-well-centered meshes. 



ial. We compared optimization of Eiq to the Mesquite implementation of Laplacian 
smoothing, applying Laplacian smoothing to the initial mesh and running it until it 
converged after 60 iterations (0.14 seconds). The result of Laplacian smoothing is a 
mesh in which 22.33% of the tetrahedra are not 3-well-centered. Figure 7.12 includes 
the h{v, a) /R{<j) distributions for the initial mesh, the mesh after optimizing i?i6, and 
the mesh resulting from Laplacian smoothing. Near each histogram we show the per- 
centage p and number n of tetrahedra (not h/R values) that are not 3-well-centered, 
and we report the mean fi and standard deviation a of the distribution of h/R values. 

It is worth noting that, because of its difficulty, obtaining well-centered triangu- 
lations and/or acute triangulations of 3-dimensional objects is significant no matter 
how they are obtained. In our other work we have made use of the optimization 
techniques developed in this paper to construct well-centered triangulations of several 
simple three-dimensional shapes [39] and to constructively prove the existence of an 
acute triangulation of the 3-dimensional cube [42] , solving an open problem mentioned 
in [19] and [7]. 

8. Conclusions and Research Questions. This paper shows that an n-well- 
centered simplex can be characterized in terms of the equatorial balls of its facets 
and uses this alternate characterization to prove that an n-well-centered mesh in M" 
is a Delaunay mesh. The paper introduces the related cost functions Eoo and Ep 
that quantify the well-centeredness of triangulations in any dimension, extending the 
function introduced in [40]. Some properties of the cost function are discussed, and 
it is shown that a cost function quantifying well-centeredness must be nonconvex. 

After introducing the cost function, the paper shows that the minmax angle tri- 
angulation is the optimal triangulation with respect to the Eoo energy and discusses 
why our algorithm uses the local preprocessing algorithm of [40] instead of comput- 
ing the maxmin triangulation after each step of optimization. The discussion raises 
the interesting research question of how to efficiently compute (and recompute) a tri- 
angulation that minimizes the maximum angle among triangulations with no lonely 
vertices. 
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The task of developing a local preprocessing algorithm that works in dimensions 
higher than 2 is another important research objective. A simple and complete char- 
acterization of the mesh connectivity requirements for a vertex and its one-ring in 
a tetrahedral mesh in M.^ to be 3-well-centered would be helpful. We have made a 
start for suc;li a characterization in [41], where we have discovered some beautiful 
connections to the triangulation of the spherical link of the one ring. 

The experiments of Section 7 show that the proposed cost function can be effective 
in finding a well-centered triangulation for meshes that permit such triangulations. 
The optimization problem in the context of our nonconvex cost functions Ep is a 
difficult problem, though, and Mesquite does not always find a global minimum of the 
energy. While it is easy to show that our gradient descent type algorithm converges to 
a local stationary point, it would be nice to have an optimization method guaranteed 
to find a global minimum of the energy. This however is a very hard problem and 
typical of the diflicnilties facied by other iterative algorithms for mesh optimization. For 
example, for the vastly popular iterative algorithms for centroidal Voronoi tessellations 
[12] and their variations [13, 14], restricted convergence results have only recently 
started appearing [11, 18]. Similarly, a convergenc;e proof for variational tetrahedral 
meshing [2] is known for only one rings, although the algorithm is very useful in 
practice. 

It would also be worthwhile to improve the efficiency of our optimization. In 

particular, it would be interesting to study methods for localizing the energy and 
applying optimization in only those specific areas where it is needed. Besides possibly 
making the optimization more efficient, localizing the energy would make it easier to 
parallelize the algorithm. The experiment in Section 7.3 that made the optimization 
easier by repositioning the boundary vertices suggests that using a constrained op- 
timization with boundary vertices free to move along the boundary could make the 
optimization more effective. 

It is also possible that the cost function could be improved. Using a linear com- 
bination of Simr with Ep was effective for the two holes mesh of Section 7.3 and the 
geographical meshes in Section 7.6, but the coefficients of the linear combination were 
chosen quite arbitrarily, and there may be other, better ways to prevent element in- 
version. There were also some experiments which needed to use Ep with more than 
one parameter p in order to find a nice result. Taking a linear combination of Ep 
for different powers of p might be effective for those situations and perhaps more 
generally. 

Since the original submission of this manuscript, the authors have become aware 
that Sazanov et al. generated a 3-well-centered mesh of a spherical layer by repeating 
the near-boundary triangulation of their mesh stitching approach without stitching to 
an ideal mesh [34] . Generalizing their construction to more complicated 3-D domains 
is another interesting direction for research. 

To summarize the paper briefly, our generalized characterization of well-centered- 
ness offers, for the first time, a direction in which planar acute triangulations may be 
generalized. More complex three dimensional experiments will have to await a better 
preprocessing and better mathematical understanding of the topological obstructions 
to well-centeredness. 

We believe we have shown enoTigh evidence in this and related publications that 
one can produce simple three dimensional well-centered tetrahedral meshes. In pla- 
nar domains, it is already possible to produce well-centered triangulations with or 
without holes and gradations, for complex domains. It is also possible to improve 
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triangulations that are already acute. Like many other successful mesh optimization 
algorithms, a convergence theory for well-centered meshing will be discovered eventu- 
ally, we hope, either by us or by other researchers. For further developments, we felt 
the need to make available the evidence that well-centered meshes are now possible 
for experiments, and that there is a useful characterization theory for such meshes. 
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